function Plot_Stationary_Dist_NPz(parms,prob_Nxz)

    phi_of_x = @(x)1./(1 + exp(-x));
    [N_mesh,P_mesh] = meshgrid(parms.grid_N(:),phi_of_x(parms.grid_x(:)));
    
    m = 2; n = 6;
    
    figure
    
    subplot(m,n,1:3)
    mesh(N_mesh,P_mesh,prob_Nxz(:,:,1)')
    xlabel('N')
    ylabel('\phi')
    title('prob(N,\phi,z), z=1')
    
    subplot(m,n,4:6)
    mesh(N_mesh,P_mesh,prob_Nxz(:,:,2)')
    xlabel('N')
    ylabel('\phi')
    title('prob(N,\phi,z), z=2')
    
    subplot(m,n,7:8)
    hold on
    plot(parms.grid_N(:),sum(sum(prob_Nxz,3),2))
    temp = squeeze(sum(prob_Nxz,2));
    plot(parms.grid_N(:),temp(:,1)/sum(temp(:,1)))
    plot(parms.grid_N(:),temp(:,2)/sum(temp(:,2)))
    xlabel('N')
    legend({'uncond','cond: z=1','cond: z=2'})
    legend boxoff
    
    subplot(m,n,9:10)
    hold on
    temp = squeeze(sum(prob_Nxz,1));
    grid_phi = phi_of_x(parms.grid_x(:));
    plot(grid_phi(:),sum(sum(prob_Nxz,3),1))
    plot(grid_phi(:),temp(:,1)/sum(temp(:,1)))
    plot(grid_phi(:),temp(:,2)/sum(temp(:,2)))
    xlabel('\phi')
    
    subplot(m,n,11:12)
    hold on
    bar(1:parms.Nz,sum(squeeze(sum(prob_Nxz,1)),1))
    xlabel('state')
    
    % cdf of lambda
    figure
    hold on
    temp = squeeze(sum(prob_Nxz,1));
    plot(grid_phi(:), cumsum(sum(sum(prob_Nxz,3),1)));
    plot(grid_phi(:),cumsum(temp(:,1)/sum(temp(:,1))))
    plot(grid_phi(:),cumsum(temp(:,2)/sum(temp(:,2))))
    grid on
    xlabel('\phi')
    title('CDF(\phi)')
    legend({'uncond','cond: z=1','cond: z=2'})
    legend boxoff
    set(gca,'ylim',[0 1])

end